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ABSTRACT 

The relationship between the boundaries for Hill and Lagrange stability in or- 
bital element space is modified in the case of resonantly interacting planets. Hill 
stability requires the ordering of the planets to remain constant while Lagrange 
stability also requires all planets to remain bound to the central star. The Hill 
stability boundary is defined analytically, but no equations exist to define the 
Lagrange boundary, so we perform numerical experiments to estimate the loca- 
tion of this boundary. To explore the effect of resonances, we consider orbital 
element space near the conditions in the HD 82943 and 55 Cnc systems. Previous 
studies have shown that, for non-resonant systems, the two stability boundaries 
are nearly coincident. However the Hill stability formula are not applicable to 
resonant systems, and our investigation shows how the two boundaries diverge 
in the presence of a mean-motion resonance, while confirming that the Hill and 
Lagrange boundaries are similar otherwise. In resonance the region of stability 
is larger than the domain defined by the analytic formula for Hill stability. We 
find that nearly all known resonant interactions currently lie in this extra stable 
region, i.e. where the orbits would be unstable according to the non-resonant Hill 
stability formula. This result bears on the dynamical packing of planetary sys- 
tems, showing how quantifying planetary systems' dynamical interactions (such 
as proximity to the Hill-stability boundary) provides new constraints on planet 
formation models. 

Subject headings: methods: N-body simulations, stars: planetary systems 



1. Introduction 

By the end of 2006, 20 multiple planetary systems had been detected beyond the Solar 
System (Butler et al. 2006, Wright et al. 2007). Of these, 7 are likely to contain at least 1 
pair that is in a mean-motion resonance. Barnes & Quinn (2004; hereafter BQ) showed that 
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one of these resonant pairs, HD 82943 b and c, had best-fit orbital elements that placed the 
system near a stability limit. Indeed, the best fit was unstable, but a small change (well 
within observational uncertainties) in the eccentricity e of the outer planet would make the 
system stable (BQ; Ferraz-Mello et al. 2005). BQ also showed that stability requires the ratio 
of the orbital periods, P c /Pb, be near 2, and that the relative mean longitudes and difference 
in longitudes of pericenter lie in a range such that conjunctions never occur at the minimum 
distance between the orbits. This result suggests that, given the values of e and a of the two 
planets, stability is only possible if the two planets are in the 2:1 resonance. 

Two types of stability have been considered in the literature. Hill stability requires 
the ordering of planets to remain constant for all time; the outer planet may escape to 
infinity. The equations that define the limits of Hill stability (i.e. Marchal & Bozis 1982; 
Gladman 1993) only apply to systems of 2 planets that are not in a resonance. Lagrange 
stability requires all planets remain bound to the star, and the orbits evolve at least quasi- 
periodically. Lagrange stability is more meaningful, but its criteria have not been delineated 
analytically. 

Barnes & Greenberg (2006a; hereafter BG) showed that the Hill-stability boundary is 
nearly the same as the Lagrange-stability boundary, at least for the non-resonant planets in 
HD 12661 and 47 UMa. Although the Hill stability boundary was derived for non-resonant 
systems, it is not clear how mean- motion resonances distort it. Here we explore the stability 
boundary near two resonant systems, HD 82943 (Mayor et al. 2004) and 55 Cnc (Marcy et 
al. 2002; McArthur et al. 2004). Note that for both systems the inner planet of the resonant 
pair is named b and the outer c. We find that the resonances do provide extra regions 
of Lagrange stability in phase space that extend beyond the analytic criterion. In § 2 we 
describe Hill and Lagrange stability and our numerical methods. In § 3 we present our 
results for HD 82943 and 55 Cnc. We also tabulate proximities to the Hill boundary for 
all applicable systems and find that all but one resonantly interacting pair would lie in an 
unstable region if not for the resonance. In § 4 we draw conclusions and suggest directions 
for future work. 



2. Methodology 



2.1. Stability Boundaries 



Hill stability in a coplanar system can be described by the following inequality: 
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where M is the total mass of the system, mi is the mass of the more massive planet, is the 
mass of the less massive planet, is the mass of the star, G is the gravitational constant, 
M* = mirri2 + m\m% + m2fn 3 , c is the total angular momentum of the system, and h is the 
energy (Marchal & Bozis 1982). If a given three-body system satisfies the inequality in Eq. 
([1]), then the system is Hill stable. If this inequality is not satisfied, then the system may or 
may not be Hill stable. In this inequality, the left-hand side is a function of the orbits, but the 
right-hand side is only a function of the masses. This approach is fundamentally different 
from other common techniques for determining stability which exploit resonance overlaps 
(Wisdom 1982; Quillen & Faber 2006), chaotic diffusion (Laskar 1990; Pepe et al. 2007), fast 
Lyapunov indicators (Froeschle et al. 1997; Sandor et al. 2007), or periodic orbits (Voyatzis 
& Hadjidemetriou 2005, 2006; Hadjidemetriou 2006). 

BG use (3 (the left-hand side of Eq. [I]), and f3 crit (the right-hand side) to describe the 
Hill stability boundary. The Hill stability boundary is the curve defined by (3/ Peru — 1- BG 
showed that the Lagrange stability boundary appears to be located where (3//3 C rit is slightly 
larger than 1 (1.02 for 47 UMa and 1.1 for HD 12661). 



2.2. Numerical Methods 

For each system, HD 82943 and 55 Cnc, 1000 initial configurations were generated based 
on the observational data for each system (Mayor et al. 2004; Marcy et al. 2002), that is, the 
initial conditions spanned the range of observational uncertainty. Note that more recent, 
improved elements are available (Butler et al. 2006), but for our purposes the older values 
serve equally well. Orbital parameters that have known errors, such as e and the period, 
P, are varied as a Gaussian centered on the best fit value, with a standard deviation equal 
to the published uncertainty, and orbital elements are sampled appropriately. For elements 
with systematic errors, such as inclination, the initial conditions were varied uniformly. 
The inclination was varied between and 5°, and the longitude of ascending node was 
varied between and 2tt. Masses were then set to the observed mass divided by the sine of 
the inclination. The variation of orbits out of the fundamental plane will not significantly 
affect our calculations of Hill stability (Veras & Armitage 2004). Each element was varied 
independently. The distribution of initial conditions is presented in Table 1. In this table, 
w is the longitude of periastron and T peri is the time of periastron passage. The integrations 
were performed with SWIFT (Levison & Duncan 1994) and MERCURY (Chambers 1999), 
and conserve energy to at least 1 part in 10 4 . For more details on these methods, consult 
BQ. 

For each simulation we numerically determine Lagrange stability on ~ 10 6 year timescales. 
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BQ showed that this timescale identifies nearly all unstable configurations. We then calcu- 
late /3/Pcrit in the parameter space sampled by the numerical integrations. Comparison of 
these two sets of results shows how the Hill and Lagrange stability boundaries are related 
near a mean-motion resonance. 



3. Results 

For HD 82943 the "stability map" is shown by the grayscale shading in Fig. [TJ Shad- 
ing indicates the fraction of initial conditions, in a certain range of orbital element space, 
that give Lagrange stable behavior (no ejections or exchanges) over 10 6 years: White bins 
contained only stable configurations, black only unstable, and the darkening shades of gray 
correspond to decreasing fractions of stable configurations. This representation plots the 
stable fraction as a function of two parameters: the eccentricity of the outer planer, e c and 
the ratio of the periods P c /Pb- The numerical simulations show that Lagrange stability is 
most likely for values of P c /Pb slightly greater than 2, and e c less than 0.4. BQ called this 
feature the "stability peninsula" . 

Superimposed on this grayscale map are contours of /3//3 C rit values. All the values of 
PI Pent are well less than 1.02 which ordinarily would imply instability. However, in the 
resonance zone where P c /Pb ~ 2, the stability peninsula sticks into a regime (/?/ /3 C rit as 
small as 0.75) that would be unstable if the planets were not in a mean motion resonance. 
Note that the numerical simulations include cases with variations of a few per cent in mass; 
the P/Pcrit contours shown are for the average mass, but would shift only slightly over the 
range of masses. 

For 55 Cnc, the stability map (Fig. [2]) was developed from integrations over 4 million 
years, i.e. 10 6 orbits of the outer planet. Our simulations include planets b, c and d, but not 
the inner planet, e. Planet e is relatively small and Zhou et al. (2004) found that the outer, 
non-resonant planet d does not appear to affect the global stability of the system. There- 
fore our simulations should elucidate the relationship between Hill and Lagrange stability 
boundaries in the presence of a 3:1 mean motion resonance. 
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Fig. 1. — Stability map for HD 82943. The shading represents the fraction of simulations that are 
Lagrange stable on 10 6 year timescales. The bin sizes are 0.02 for e c and 0.01 for P c /Pb- Contour 
lines show values of f3/f3 cr it- Ordinarily (3//3 cr it < 1-02 would imply instability, however the mean 
motion resonance provides a stable region that would not exist if the resonance did not affect the 
motion. 
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Fig. 2. — Stability map for 55 Cnc (c.f. Fig. [[]). The bin sizes are 0.05 for e c and 0.04 for P c /Pb- 
As in HD 82943, the mean-motion resonance provides a larger Lagrange-stable region. 
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Of our simulations 50.2 ± 5.5% were Lagrange stable. The least massive planet, c, was 
the planet most likely to be ejected. In this system we see that stability is likely for e c < 0.3 
everywhere, except at P c /Pb ~ 3, where it extends to e ~ 0.55. 

Comparing this distribution with the analytical stability criterion (/3//3 C rit <l-03) we see 
that the numerical experiments reproduce that boundary, except in resonance where (3 / /3 cr it 
can be as low as 0.96. This stability peninsula for 55 Cnc does not protrude as far into 
the Hill unstable region as HD 82943. This difference may be because the 2:1 resonance 
is of a lower order (and thus stronger) than the 3:1, and therefore has a more pronounced 
stabilizing effect. 

Next we tabulate (3/(3cHt values for all observed systems that contain two planets. We 
also include GJ 876 c and b, and v And c and d. Eq. (jTJ is only applicable to two-planet 
systems, but we consider these latter two pairs, which are each part of a bigger system, as 
the third planet in each system is probably too small or too far away to significantly change 
the interaction of those pairs. However, interpreting values of (3 / ' (3 CT it in systems of more 
than two companions should be made with caution, as there is no guarantee f3/(3 cr n — 1 
corresponds to the Hill boundary for any individual pair. 

Table 2 lists values of (3/f3 cr u and the "class" of the interaction, which distinguishes 
the dominant phenomenon that changes the shapes of the orbits. "R" denotes pairs whose 
interaction is dominated by mean-motion resonances (Table 2 also lists the resonance), "T" 
indicates pairs that may have experienced significant tidal evolution, and "S" indicates pairs 
with strong secular interactions (Barnes & Greenberg 2006b). All but one resonantly inter- 
acting pair have (3/ ' j3 CT it values less than 1. If not for the resonance, these systems would be 
unstable. 

Overall, we find 70% of the pairs we consider are observed to have (3 / (3 crit < 1.3. HD 
217107 is observed to have a (3/ f3 C rit significantly larger than other pairs. In Fig. [3] we plot 
the distribution of /3//3 C rit values from Table 2. 



- 8- 



0.25 



"i i i i i i i i i i i i i i r 



0.2 



0.15 



o 
. , — i 
-+-> 

o 



0.1 



0.05 







J I L 



J L 



2 



J L 



4 

ft/Per* 



J L 



6 



8 



Fig. 3. — Histogram of values of (3/(3 C rit (with a bin size of 0.1) for 17 pairs of planets in Table 2. 
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4. Conclusions 

By explicitly mapping how mean-motion resonances can provide additional regions of 
stability in orbital element space, we have found that nearly all observed resonant systems 
lie in these extended regions. More generally, we have also shown the distribution of (3/ f3 cr i t 
appears to show that many planetary systems (resonant or not) lie close to the limits of 
dynamical stability. These distributions provide new constraints for models of planet forma- 
tion. 

In the cases presented here, the 2:1 resonance provides a larger stable region than 3:1, 
presumably because it is a lower order (stronger) resonance. However, for the 5:1 mean- 
motion resonance in HD 202206, (3/ 'Peru can reach 0.88 and still be stable. So why is the 
range of stability for the 3:1 resonance in 55 Cnc so small? Perhaps if e c in the 55 Cnc 
system has values in excess of 0.5, it does interact with the third planet, destabilizing the 
system. Future work should investigate the minimum P/P cr it that allows stability for each 
resonance. Future work may also determine if the peninsula we find in 55 Cnc is truncated 
due to interactions with 55 Cnc d. 

We seek to identify the origin of the shape of the stability peninsulas in resonant systems. 
Ideally a general expression will eventually be developed that describes the Lagrange stability 
boundary that applies to planets both in and out of resonance. One avenue of research is to 
focus on close approach distances. In the limit of zero eccentricity, orbits are unstable if they 
are separated by less than 3.5 mutual Hill radii (Gladman 1993). Therefore we speculate that 
systems with approaches within this distance are unstable. For secularly evolving systems, 
the orbits change with time and eventually the planets will line up at the minimum distance 
between the orbits over a secular period. Resonances can prevent planets from lining up at 
this danger zone, hence the stability peninsulas. This likely explanation for the shape of the 
Lagrange boundaries might be a fruitful direction of future research into planetary system 
stability. 

The distribution of P/P cr it shows that, regardless of the the presence of mean-motion 
resonance, many systems have values that are close to the stability boundary. This trend 
appears to support the hypothesis that planetary systems are dynamically "packed", i.e. 
that additional planets could not exist in orbits between those that are known without 
destabilizing the system (Barnes & Quinn 2001; BQ, Barnes & Raymond 2004; Raymond & 
Barnes 2005; Raymond et al. 2006). Perhaps there is a minimum value of P/Pcrit that would 
permit the insertion of an additional planet that leaves the system still stable. In other 
words it will be interesting to determine, for a given value of P/Pcru, the largest mass object 
that could orbit between two planets and still leave the system stable. Such a relation could 
produce an analytic criterion for dynamical packing, which can currently only be estimated 
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numerically (e.g. Menou & Tabachnik 2003; Dvorak et al. 2003; Rivera & Haghighipour 
2007). 

Past work provides illumination on the possibility that some minimum value of (5/ f3 crit 
may define the limit for which additional planets could be placed between the observed 
planets. The HD 168443 system (P/f3 C rit = 1-94) has been shown to be unable to support 
even infinitesimal test masses (Barnes & Raymond 2004). The region between the known 
planets of HD 169830 (fi/ficrit = 1-28) is chaotic and a planet in that region is most likely 
unstable (Erdi et al. 2004). On the other hand, HD 38529 ((3/(3 crit = 2.07) could support 
a Saturn- mass companion between the known planets (Barnes & Raymond 2004). These 
results suggest fi/ficrit = 2 may be the critical value. 

The packing of the two planets in HD 190360 demands closer inspection. Although the 
orbits are more separated and less eccentric than those in HD 168443, their f3/f3 cr i t value 
(1.70) is less than that for HD 168443 (1.94). To explore this issue, we have integrated the 
HD 190360 system with a hypothetical Earth-mass planet on a circular orbit located at the 
midpoint between the apoastron distance of the inner planet and the periastron distance 
of the outer. The additional companion in the HD 190360 system survived for 10 6 years. 
A similar experiment with HD 74156 (f3/f3 crit = 1.54) showed ejection of the Earth-mass 
planet in only 2500 years. We tentatively conclude that systems are packed if fi/ficrit ^1-5, 
not packed if f3/f3 cr it >2, and the packing status is unknown in the range 1.5 <f3/f3 C rit <2. 
Future research needs to determine the relationship between (3/(3 C rit and the possibility that 
additional planets could be stable between known ones. 

As noted at the end of § 3.3, 70% of the tabulated systems have f3 / f3 crit < 1.3, indi- 
cating that the planets are too fully packed to allow any intermediate planets. This result, 
coupled with the limitations of radial velocity surveys to detect planets (e.g. we used mini- 
mum masses), suggests that the vast majority of multiple planet systems are similarly fully 
packed. Our results are therefore consistent with the "Packed Planetary Systems" hypothe- 
sis (BQ; Barnes & Raymond 2004; Raymond & Barnes 2005; Raymond et al. 2006; see also 
Laskar 1996) which proposes that planetary systems tend to form so as to be dynamically 
packed. This hypothesis therefore predicts that HD 190360 and especially HD 217107 harbor 
additional, undetected planets. 

This investigation has identified a simple way to parameterize multiple planet systems. 
At least for a two-planet system, a single parameter {3//3 cr it may indicate both stability 
and packing. Moreover, the statistics of the distribution of this dynamical parameter for 
observed systems are intriguing: Planetary systems tend to be dynamically fully packed 
and resonant systems lie at values of (3/ f3 cr u that would indicate instability for non-resonant 
systems. Describing planetary systems by parameterizing the character of their dynamical 
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interaction is also the approach taken by Barnes & Greenberg (2006b), who calculated the 
proximities of planetary systems to the apsidal separatrix. These new methodologies focus 
on the proximities of the dynamical interactions to boundaries between qualitatively different 
dynamical regimes. 

It now appears that about half of stars with planets have multiple planets (Wright et 
al. 2007), and descriptions of their dynamical interactions will therefore become increasingly 
more relevant, especially since many planets' eccentricities oscillate by two orders of mag- 
nitude (Barnes & Greenberg 2006b). We encourage research that models planet formation 
(e.g. Lee & Peale 2002; Sandor & Kley 2006) to include comparisons of the simulated values 
of fl/Pcrit to those of real planetary systems. 

J. Bryan Henderson, Thomas R. Quinn, and Chance Reschke assisted with the simu- 
lations presented here. An anonymous referee provided helpful suggestions. This work was 
funded by NASA's PG&G program. 
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Table 1: Orbital Elements and Errors for HP 82943 and 55 Cnc 



System m 3 (M ) Planet m (M Jup ) 



P(d) 



- pen 



(JD) 



HD 82943 1.05±0.05 a 



55 Cnc 0.95±0.1 b 



b 

c 
b 

c 

d 



0.88 

1.63 
0.84 ±0.07 
0.21 ±0.04 
4.05 ±0.9 



221.6 ±2.7 
444.6 ± 8.8 
14.653±0.0006 
44.276±0.021 
5360±400 



0.54 ±0.05 
0.41 ±0.08 
0.02±0.02 
0.339±0.21 
0.16±0.06 



138 ± 13 
96 ±7 
99±35 
61±25 
201±22 



2451630.9 ± I 
2451620.3 ± 1 
2450001.479±: 
2450031.4±2 
2785±250 



a Santos et al (2000); b Marcy et al. (2002) 



Table 2: Values of (3/ (3 crit for Known Systems 
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System 


Pair 




Class 


HD 202206° 


b-c 


0.883 


R (5:1) 


HD 128311 6 


b-c 


0.968 


R (2:1) 


HD 82943 b 


b-c 


0.946 


R(2:l) 


HD 73526 6 


b-c 


0.982 


R(2:l) 


GJ 876 6 


c-d 


0.99 c 


R(2:l) 


47 UMa 6 


b-c 


1.025 


S 


HD 155358 d 


b-c 


1.043 


S 


HD 108874 6 


b-c 


1.107 


R (4:1) 


v And 6 


c-d 


1.125 c 


S 


HD 12661 6 


b-c 


1.199 


s 


HIP 14810 6 


b-c 


1.202 


T 


HD 169830 6 


b-c 


1.280 


S 


HD 74156 6 


b-c 


1.542 


s 


HD 190360 6 


b-c 


1.701 


T 


HD 168443 6 


b-c 


1.939 


S 


HD 38529 b 


b-c 


2.070 


s 


HD 217107 6 


b-c 


7.191 


T 



° Correia et al. (2005); b Butler et at. (2006); c An additional planet is present in this system; 
d Cochran et al. (2007); e Wright et al. (2007) 



